

subroutine  compute_Richter_ML(confidence,c,nd,stanum,evnum,nsta,nev,sigma2_RH)
    implicit none
    
    include 'seisev_EVENT_PHASE_structure.FI'
    include 'inversione_pesata_multipar.FI'
    
    integer i,k,nd,isw,nv,iter,n_par,lun
    integer stanum(150000),evnum(150000),nsta,nev
    integer RichDist(100),first,last,dst
    real*4  RichVal(100),dif,corr,resmax,confidence
    real*4  v(5000),wmag,wsd,wmax,wmin,ww(5000),ruse
    real*8  sigma2_RH
    real*4  stares_vet(5000),sta_count_vet(5000),rms
    real*8  sig2_old,sig2_new,qtot,signif
    logical significant
    character*5 memosta

    data isw/0/

   !
    n_par = nsta + nev
    
    do i = 1,nd
      c(i).MLRH_corr = 0.0
    enddo
  ! compute (not read) epicentral distances from hypocentral and depth.
    do i=1,nd
       c(i).MLRH_dist = SQRT(  c(i).s_h_hypdist * c(i).s_h_hypdist -  &
                               c(i).depth       * c(i).depth )
    enddo


   ! compute all the station Richter magnitudes
    if(isw .eq. 0)then
       isw = 1
       call read_Richter_tab(nv,RichDist,RichVal)
    endif

    iter=0
    significant = .TRUE.
1   do while(significant .and. iter .lt. 100 )
    iter = iter +1
    
    do i = 1,nsta
      stares_vet(i)=0.0
      sta_count_vet(i) = 0.0
    enddo
 
    
    do i=1,nd       
       dst = NINT(c(i).MLRH_dist)
       ! compute correction
       call compute_Richter_corr(dst,corr,nv,RichDist,RichVal)
       ! compute station magnitude
       c(i).MLRH_sta = log10(c(i).WA_pha_amp) + corr + c(i).MLRH_corr 
!       write(*,*)evnum(i),c(i).sta,c(i).MLRH_sta,c(i).MLDB_sta,c(i).MLDB_ev

    enddo
    
!   Prima ricalcolo le magnitudo di evento con la mediana,
!   poi con la media di Huber (tanto per vedere la differenza
!   delle due varianze).


    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1

    sigma2_RH=0.0d00
    do while (last .le. nd  .and. first .le. nd)

      k= 0
      do i=first,last
          k=k+1 
          v(k) = c(i).MLRH_sta
      enddo

      call mediana(v,k,wmag)      

      k=0
      do i=first,last
          k=k+1
          c(i).MLRH_ev = wmag     !event mag
          dif = c(i).MLRH_sta - wmag
          sigma2_RH = sigma2_RH + dif*dif
      enddo   
                 
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1

    enddo        
 
    write(*,*)' calcolo delle ML Richter con la mediana. Sigma2 = ',sigma2_RH    


    ! go throught the events and compute the event Richter Magnitudes.
    ! ora calcolo le magnitudo RIchter finali con Huber

    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1

    sigma2_RH=0.0d00
    do while (last .le. nd  .and. first .le. nd)

      k= 0
      do i=first,last
          k=k+1 
          v(k) = c(i).MLRH_sta
      enddo

      ruse = 0.23
      wmag = 0.
      wsd = 0.
      call Huber_dev(k,v,ruse,wmag,wsd,wmax,wmin,ww)
      ! memorizza i risultati   !
      !
      
      k=0
      do i=first,last
          k=k+1
          c(i).MLRH_ev = wmag     !event mag
          dif = c(i).MLRH_sta - wmag
          c(i).MLRH_s_wgt = ww(k)           ! memorizzo qui i pesi Richter

          c(i).MLRH_dif = dif

          sigma2_RH = sigma2_RH + dif*dif
!
          stares_vet(stanum(i)) = stares_vet(stanum(i)) + dif*ww(k) 
          sta_count_vet(stanum(i)) = sta_count_vet(stanum(i)) +ww(k) 
      enddo   
    
    
                 
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1

    enddo        


    write(*,*)' calcolo delle ML Richter con Huber. Sigma2 = ',sigma2_RH  
    
      
      rms=SNGL(sigma2_RH)
      rms=rms/float(nd)
      rms=SQRT(rms)
       write(11,'(a,a,i3,a,f10.2,a,i5,a,f8.4)') &
       'Richter ',',',iter,',',sigma2_RH,',',nd,',',rms
      
        
    ! i residui secondo RH
    do i = 1,nsta
      stares_vet(i) = stares_vet(i) / sta_count_vet(i)
    enddo
    
    resmax = -1.
    do i = 1,nd
       c(i).MLRH_corr = c(i).MLRH_corr - stares_vet(stanum(i))
       if(ABS(stares_vet(stanum(i))) .gt. resmax)then
         resmax = stares_vet(stanum(i))
         memosta = c(i).sta
       endif
    enddo 
    write(*,*)' residuo ABS massimo ',memosta,' ',resmax 
    if(iter.eq.1)then
       sig2_old = sigma2_RH
    else
       sig2_new = sigma2_RH
       call Ftest(nd,sig2_old,n_par,sig2_new,n_par,qtot,signif)
       sig2_old = sigma2_RH
       if(signif .gt. confidence  .and. resmax .gt. 0.005)then
          significant = .TRUE.
       else
          significant = .FALSE.
       endif
       write(*,*)' Richter residuals; iter signif resmax',iter,signif,resmax
       
    endif
    enddo 
    
! calcolo finale della variazione dei valori di ML di stazione per ogni stazione

    do i=1,nsta
      stares_vet(i) = 0.
      sta_count_vet(i) = 0.
    enddo
    do i=1,nd
      stares_vet(stanum(i)) = stares_vet(stanum(i)) +  &
                              c(i).MLRH_dif*c(i).MLRH_dif
      sta_count_vet(stanum(i)) = sta_count_vet(stanum(i)) + 1.
    enddo
    do i=1,nsta
      stares_vet(i) = SQRT( stares_vet(i) / (sta_count_vet(i)-1.))
    enddo
    
! OUTPUT DEI RISULTATI RICHTER ============================================
! *************************************************************************
! =========================================================================

    call get_lun(lun)
    open(unit=lun,file='results_RICHTER.csv',status='unknown')
    write(lun,'(a)')'Richter,iter,sigma2_HB,nd,rms' 
       write(lun,'(a,a,i3,a,f10.2,a,i5,a,f8.4)') &
       'Richter ',',',iter,',',sigma2_RH,',',nd,',',rms
    


    write(lun,'(a)')'Sta-code,corr,sigma'
    do i=1,nsta
      do k=1,nd
        if(stanum(k).eq.i)then
          write(lun,'(a,2(a,f7.4))')c(k).sta,',',c(k).MLRH_corr,',',stares_vet(i)
          goto 99
        endif
      enddo
 99   continue
    enddo
    close(lun)    
return
end

    
 